Abstract
Include your Abstract here. This should be an unstructured abstract of 200 words, summarising the main contents of the scientific report.This document provides the code used to carry out the analysis of the project on Bayesian modelling of crime against women in India, for reproducibility purposes. The code is arranged in the following sections, which can be navigated using “Outline” in RMarkdown:
Descriptive analysis
Modelling: this section includes code used to fit three models: spatial, spatio-temporal, and spatio-temporal with Interaction type I, each fitted separately for rape and dowry deaths, as represented in the subtitles. It also includes the code for visualising the results of the models.
Sensitivity analysis
Results visualisation: the code used to produce additional figures and compile them into the visualisation file that was used in the report.
# loading all required libraries
library(dplyr)
library(sf)
library(spdep)
library(tidyr)
library(INLA)
library(ggplot2)
library(viridis)
library(patchwork)
library(knitr)
library(kableExtra)
library(scales)
library(sp)
library(grid)
# Setting up the workspace
rm(list=ls())
#setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
load("CrimeUttarPradesh.RData")
carto_up <- st_as_sf(carto_up)
ls()
## [1] "carto_india" "carto_up" "data"
# Calculate crude rates per 100,000 women - to be used in dataframe below for Supplementary Figure 1
data$rape_rate_per100k <- data$rape / data$pop * 100000
data$dowry_rate_per100k <- data$dowry/ data$pop * 100000
# Create data aggregated across time (i.e. collapsed across space)
year_agg <- data %>% group_by(year) %>% summarise(rape_obs = sum(rape),
rape_exp = sum(e_rape),
dowry_obs = sum(dowry),
dowry_exp = sum(e_dowry),
total_pop = sum(pop),
SMR_rape = sum(rape)/sum(e_rape),
SMR_dowry = sum(dowry)/sum(e_dowry),
rape_crude_rate = rape_obs / total_pop * 100000,
dowry_crude_rate = dowry_obs / total_pop * 100000,
max_rape_region = dist[which.max(rape_rate_per100k)],
max_rape_rate = max(rape_rate_per100k),
min_rape_region = dist[which.min(rape_rate_per100k)],
min_rape_rate = min(rape_rate_per100k),
max_dowry_region = dist[which.max(dowry_rate_per100k)],
max_dowry_rate = max(dowry_rate_per100k),
min_dowry_region = dist[which.min(dowry_rate_per100k)],
min_dowry_rate = min(dowry_rate_per100k))
# Visualise temporal trends of SMRs
usmoothed_map <- ggplot(year_agg, aes(x = year)) +
geom_line(aes(y=SMR_rape, color = "Rape")) +
geom_line(aes(y=SMR_dowry, color = "Dowry")) + theme_minimal() +
labs(title = "Temporal evolution of SMRs of Rape and Dowry Crime",
x = "Year", y = "Crude Rate", color = "Crime Type")
usmoothed_map
# Visualise summary stats - regions with max & min crude rates each year
year_agg_vis <- year_agg %>% mutate(
'Year' = year,
'Rape' = round(rape_crude_rate, 1),
'Dowry' = round(dowry_crude_rate, 1),
'Higest Rape' = paste0(max_rape_region, " (", round(max_rape_rate, 1), ")"),
'Lowest Rape' = paste0(min_rape_region, " (", round(min_rape_rate, 1), ")"),
'Highest Dowry'= paste0(max_dowry_region, " (", round(max_dowry_rate, 1), ")"),
'Lowest Dowry' = paste0(min_rape_region, " (", round(min_dowry_rate, 1), ")") ) %>%
select( 'Year', 'Rape', 'Dowry', 'Higest Rape', 'Lowest Rape', 'Highest Dowry', 'Lowest Dowry')
reg_table <- knitr::kable(year_agg_vis, caption = "Yearly summary of crude crime rates (per 100,000 women)") %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center")
reg_table
| Year | Rape | Dowry | Higest Rape | Lowest Rape | Highest Dowry | Lowest Dowry |
|---|---|---|---|---|---|---|
| 2001 | 5.4 | 6.1 | Pilibhit (18.4) | Sant Ravidas Nagar Bhadohi (0.3) | Kheri (12.9) | Sant Ravidas Nagar Bhadohi (1.5) |
| 2002 | 3.8 | 5.1 | Pilibhit (11.3) | Sant Ravidas Nagar Bhadohi (0) | Etah (14.7) | Sant Ravidas Nagar Bhadohi (1.2) |
| 2003 | 2.4 | 3.4 | Hathras (7.2) | Sant Kabir Nagar (0) | Mainpuri (8.7) | Sant Kabir Nagar (0.8) |
| 2004 | 3.5 | 4.3 | Sitapur (8.6) | Ghazipur (0.6) | Firozabad (10) | Ghazipur (0.8) |
| 2005 | 3.0 | 3.8 | Aligarh (8.4) | Ghazipur (0.3) | Firozabad (9.7) | Ghazipur (0.4) |
| 2006 | 3.1 | 4.3 | Gautam Buddha Nagar (7.7) | Sant Ravidas Nagar Bhadohi (0.6) | Kanpur Dehat (13.4) | Sant Ravidas Nagar Bhadohi (1.6) |
| 2007 | 3.8 | 4.8 | Chitrakoot (12.1) | Deoria (0.3) | Mainpuri (12.2) | Deoria (1.5) |
| 2008 | 4.2 | 5.1 | Chitrakoot (11.2) | Sant Ravidas Nagar Bhadohi (0.6) | Aligarh (11.1) | Sant Ravidas Nagar Bhadohi (1.1) |
| 2009 | 3.9 | 4.9 | Chitrakoot (9.9) | Ghazipur (0.4) | Aligarh (10.2) | Ghazipur (1.9) |
| 2010 | 3.3 | 4.8 | Rampur (12) | Sant Ravidas Nagar Bhadohi (0.3) | Aligarh (11.3) | Sant Ravidas Nagar Bhadohi (1.1) |
| 2011 | 4.3 | 4.9 | Aligarh (10.2) | Sant Ravidas Nagar Bhadohi (0.5) | Etawah (11.2) | Sant Ravidas Nagar Bhadohi (1.3) |
| 2012 | 4.1 | 4.6 | Chitrakoot (14.2) | Mirzapur (0.7) | Aligarh (11) | Mirzapur (1) |
| 2013 | 6.2 | 4.7 | Hamirpur (19.3) | Mirzapur (0.8) | Firozabad (10.7) | Mirzapur (1.2) |
| 2014 | 6.8 | 4.9 | Banda (17.3) | Shrawasti (2.1) | Etah (11.8) | Shrawasti (1.3) |
Produced for supplementary figure 1
# Define color map
color_map <- c(
"Dowry: State Average" = "#8E44AD",
"Rape: State Average" = "#21918C",
"Individual Districts" = "lightgray"
)
# Dowry crime
p_dowry <- ggplot() +
geom_line(data=data, aes(x=year, y=dowry_rate_per100k, group=ID_area),
color = "lightgray", alpha=0.7, show.legend = FALSE) +
geom_line(data=year_agg, aes(x=year, y=dowry_crude_rate, color = "Dowry: State Average"), size = 1.5) +
scale_color_manual(name = "Legend", values = color_map) +
theme(legend.position = "right") +
labs( y = "Crude Rate of Dowry Crime", x = "Year") +
theme_minimal()
# Rape
p_rape <- ggplot() +
geom_line(data=data, aes(x=year, y=rape_rate_per100k, group=ID_area),
color = "lightgray", alpha=0.7, show.legend = FALSE) +
geom_line(data=year_agg, aes(x=year, y=rape_crude_rate, color = "Rape: State Average"), size = 1.5) +
scale_color_manual(name = "Legend", values = color_map) +
theme(legend.position = "right") +
labs(y = "Crude Rate of Rape", x = "Year") +
theme_minimal()
# Compule into single plot
region_temp_smrs <- ((p_dowry | p_rape) +
plot_layout(guides = "collect") +
plot_annotation(title = "Crude Rates of Rape and Dowry Crime in Uttar Pradesh (2001–2014)"))
region_temp_smrs
# Aggregate data across space & calculate required summary stats
dist_rates <- data %>%
group_by(ID_area) %>%
summarise(
rape_obs = sum(rape),
rape_exp = sum(e_rape),
dowry_obs = sum(dowry),
dowry_exp = sum(e_dowry),
SMR_rape = sum(rape)/sum(e_rape),
SMR_dowry = sum(dowry)/sum(e_dowry))
# Join with SP object
map_rates <- left_join(carto_up, dist_rates, by = c("ID_area" = "ID_area"))
# Convert to long format for easier plotting
map_rates_long_SMR <- map_rates %>%
select(ID_area, geometry, SMR_rape, SMR_dowry) %>%
pivot_longer(cols = c(SMR_rape, SMR_dowry), names_to = "CrimeType", values_to = "SMR") %>%
mutate(CrimeType = recode(CrimeType, rape_rate = "Rape",dowry_rate = "Dowry")
)
# Plot SMRs
unsmoothed_smr_map <- ggplot(map_rates_long_SMR) +
geom_sf(aes(fill = SMR), color = NA) +
scale_fill_viridis_c(
name = "SMR",
labels = label_comma(accuracy = 1)) +
facet_wrap(~ CrimeType) + # set ncol=1 argument if want top and bottom
labs(title = "SMRs by District (2001–2014)") + theme_bw()
unsmoothed_smr_map
# 1 - Check the coordinates & min distances between polygons:
#st_coordinates(st_centroid(carto_up)) # these appear to be latitudes & longitudes
coord <- st_coordinates(st_centroid(carto_up))
dists <- spDists(coord)
#min(dists[dists > 0]) # min distance is 0.309
# 2- define neighbours list & convert to adjacency matrix
carto_up$reid <- carto_up$ID_area
carto_up_nb <- poly2nb(carto_up, snap=0.001, queen=TRUE) # snap=0.001 seems the most accurate geographically
summary(carto_up_nb)
## Neighbour list object:
## Number of regions: 70
## Number of nonzero links: 344
## Percentage nonzero weights: 7.020408
## Average number of links: 4.914286
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9
## 2 1 12 14 16 14 6 3 2
## 2 least connected regions:
## Uttar Pradesh_Lalitpur Uttar Pradesh_Saharanpur with 1 link
## 2 most connected regions:
## Uttar Pradesh_Budaun Uttar Pradesh_Fatehpur with 9 links
nb2INLA("map_adj",carto_up_nb) # create object with the location of the graph
adj = inla.read.graph(filename="map_adj") # store graph
# create dataset collaped across time
up_agg <- data %>% group_by(ID_area) %>% # aggregate over areas
summarise(O = sum(dowry), E = sum(e_dowry))
# fit hierarchical poisson log-linear model
ID <- seq(1,70)
# define formula
f_BYM2 <- O ~ f(ID, model="bym2", graph=adj,
hyper=list(
prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
phi=list(prior="pc", param=c(0.5,2/3))))
# fit model
s_dowry_model = inla(formula=f_BYM2, family="poisson", data=up_agg, E=E,
control.predictor = list(compute=TRUE),
control.compute = list(waic=TRUE))
# obtain posterior RRs
RR_s <- c()
for(i in 1:70) {RR_s[i] <- inla.emarginal(function(x) exp(x), s_dowry_model$marginals.random$ID[[i]])}
# obtain posterior probabilities
RR_s_marg <- s_dowry_model$marginals.random$ID[1:70]
PP_s <- lapply(RR_s_marg, function(x) {1-inla.pmarginal(0,x)})
# Combine RRs & PPs into a dataframe
s_RR_PP <- data.frame(resRR=RR_s, PP=unlist(PP_s), SP_ID=up_agg[,1])
# Plot posterior RRs & PPs:
# 1 - Categorise variables
breaks_rr = c(min(s_RR_PP$resRR), 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(s_RR_PP$resRR))
s_RR_PP$resRR_cat = cut(s_RR_PP$resRR, breaks=breaks_rr, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
s_RR_PP$PP_cat = cut(s_RR_PP$PP, breaks_pp, include.lowest = TRUE)
# 2 - join with SP
map_s_RR_PP=left_join(carto_up, s_RR_PP, by=c("ID_area"="ID_area"))
# 3 - plot RRs
dowry_RR_s = ggplot() + geom_sf(data=map_s_RR_PP) + aes(fill=resRR_cat) + theme_bw() +
scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatial model") +
theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
# 4 - plot PPs
dowry_PP_s = ggplot() + geom_sf(data=map_s_RR_PP) + aes(fill=PP_cat) + theme_bw() +
scale_fill_viridis(
option = "plasma", name="PP",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP Spatial model") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
# create dataset collaped across time
up_agg_rape <- data %>% group_by(ID_area) %>% # aggregate over areas
summarise(O = sum(rape), E = sum(e_rape))
# fit hierarchical poisson log-linear model
ID <- seq(1,70)
f_BYM2 <- O ~ f(ID, model="bym2", graph=adj,
hyper=list(
prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
phi=list(prior="pc", param=c(0.5,2/3))))
s_rape_model = inla(formula=f_BYM2, family="poisson", data=up_agg_rape, E=E,
control.predictor = list(compute=TRUE),
control.compute = list(waic=TRUE))
# obtain posterior RRs
RR_s <- c()
for(i in 1:70) {RR_s[i] <- inla.emarginal(function(x) exp(x), s_rape_model$marginals.random$ID[[i]])}
# obtain posterior PPs
RR_s_marg <- s_rape_model$marginals.random$ID[1:70]
PP_s <- lapply(RR_s_marg, function(x) {1-inla.pmarginal(0,x)})
# combine into a dataframe
s_RR_PP <- data.frame(resRR=RR_s, PP=unlist(PP_s), SP_ID=up_agg[,1])
# Plot - as detailed in the chunk above
breaks_rr = c(min(s_RR_PP$resRR), 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(s_RR_PP$resRR))
s_RR_PP$resRR_cat = cut(s_RR_PP$resRR, breaks=breaks_rr, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
s_RR_PP$PP_cat = cut(s_RR_PP$PP, breaks_pp, include.lowest = TRUE)
map_s_RR_PP=left_join(carto_up, s_RR_PP, by=c("ID_area"="ID_area"))
rape_RR_s = ggplot() + geom_sf(data=map_s_RR_PP) + aes(fill=resRR_cat) + theme_bw() +
scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatial model") +
theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
rape_PP_s = ggplot() + geom_sf(data=map_s_RR_PP) + aes(fill=PP_cat) + theme_bw() +
scale_fill_viridis(
option = "plasma", name="PP",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP Spatial model") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
# Prepare data
data_st_dowry = left_join(carto_up, data, by = "ID_area")
data_st_dowry = data_st_dowry %>% dplyr::rename(O=dowry, E=e_dowry)
# Define formula & fit the model
formula_st_dowry = O ~ f(ID_area, model="bym2", graph=adj, hyper=list(
prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
phi=list(prior="pc", param=c(0.5,2/3))
)) +
f(ID_year, model="rw1", hyper=list(prec=list(prior="pc.prec", param=c(0.5/0.31, 0.01))))
st_dowry_model = inla(formula=formula_st_dowry, family="poisson", data=data_st_dowry, E=E, control.compute=list(waic=TRUE),
control.predictor = list(compute=TRUE))
# Derive residual posterior means & CIs of random walk component
# Note: they are not plotted yet as they will be first combined with ST+I resRRs - please see "Model Visualisation"
RR_stRW_RR = c()
RR_stRW_l = c()
RR_stRW_h = c()
for(i in 1:14) {
RR_stRW_RR[i] = inla.emarginal(function(x) exp(x), st_dowry_model$marginals.random$ID_year[[i]])
RR_stRW_l[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), st_dowry_model$marginals.random$ID_year[[i]]))
RR_stRW_h[i] = inla.qmarginal(0.975,inla.tmarginal(function(x) exp(x), st_dowry_model$marginals.random$ID_year[[i]]))
}
RR_stRW_dowry = data.frame(RR=RR_stRW_RR, low_CI=RR_stRW_l, up_CI = RR_stRW_h)
# Derive & plot RRs & PPs of the BYM2 component
RR_stBYM = c()
for(i in 1:70) {RR_stBYM[i] = inla.emarginal(function(x) exp(x), st_dowry_model$marginals.random$ID_area[[i]])}
RR_stBYM_marg = st_dowry_model$marginals.random$ID_area[1:70]
PP_stBYM = lapply(RR_stBYM_marg, function(x) {1-inla.pmarginal(0,x)})
resRR_PP_st = data.frame(RR=RR_stBYM, PP=unlist(PP_stBYM), ID_area=carto_up$ID_area)
# plot (using same process as for spatial plots above)
breaks = c(min(resRR_PP_st$RR), 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(resRR_PP_st$RR))
resRR_PP_st$RR_cat = cut(resRR_PP_st$RR, breaks=breaks, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
resRR_PP_st$PP_cat = cut(resRR_PP_st$PP, breaks_pp, include.lowest = TRUE)
map_st_RR_PP=left_join(carto_up, resRR_PP_st, by=c("ID_area"="ID_area"))
# Plot RRs
dowry_RR_st = ggplot() + geom_sf(data=map_st_RR_PP) + aes(fill=RR_cat) + theme_bw() +
scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatio-temporal model") +
theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
# Plot PPs
dowry_PP_st = ggplot() + geom_sf(data=map_st_RR_PP) + aes(fill=PP_cat) + theme_bw() +
scale_fill_viridis(
option = "plasma", name="PP",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP Spatio-temporal model") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
dowry_RR_st | dowry_PP_st
# Prep data
data_st_rape = left_join(carto_up, data, by = "ID_area")
data_st_rape = data_st_rape %>% dplyr::rename(O=rape, E=e_rape)
# Define & fit the model
formula_st_rape = O ~ f(ID_area, model="bym2", graph=adj, hyper=list(
prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
phi=list(prior="pc", param=c(0.5,2/3))
)) +
f(ID_year, model="rw1", hyper=list(prec=list(prior="pc.prec", param=c(0.5/0.31, 0.01))))
st_rape_model = inla(formula=formula_st_rape, family="poisson", data=data_st_rape, E=E, control.compute=list(waic=TRUE),
control.predictor = list(compute=TRUE))
# Repeat the process above for plotting of spatial RRs & PPs and (eventually) of temporal RRs
RR_stRW_RR = c()
RR_stRW_l = c()
RR_stRW_h = c()
for(i in 1:14) {
RR_stRW_RR[i] = inla.emarginal(function(x) exp(x), st_rape_model$marginals.random$ID_year[[i]])
RR_stRW_l[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), st_rape_model$marginals.random$ID_year[[i]]))
RR_stRW_h[i] = inla.qmarginal(0.975,inla.tmarginal(function(x) exp(x), st_rape_model$marginals.random$ID_year[[i]]))
}
RR_stRW_rape = data.frame(RR=RR_stRW_RR, low_CI=RR_stRW_l, up_CI = RR_stRW_h)
RR_stBYM = c()
for(i in 1:70) {RR_stBYM[i] = inla.emarginal(function(x) exp(x), st_rape_model$marginals.random$ID_area[[i]])} ##
RR_stBYM_marg = st_rape_model$marginals.random$ID_area[1:70]
PP_stBYM = lapply(RR_stBYM_marg, function(x) {1-inla.pmarginal(0,x)})
resRR_PP_st = data.frame(RR=RR_stBYM, PP=unlist(PP_stBYM), ID_area=carto_up$ID_area)
# plot (using same process as described above)
breaks = c(min(resRR_PP_st$RR), 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(resRR_PP_st$RR))
resRR_PP_st$RR_cat = cut(resRR_PP_st$RR, breaks=breaks, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
resRR_PP_st$PP_cat = cut(resRR_PP_st$PP, breaks_pp, include.lowest = TRUE)
map_st_RR_PP=left_join(carto_up, resRR_PP_st, by=c("ID_area"="ID_area"))
# RR plot
rape_RR_st = ggplot() + geom_sf(data=map_st_RR_PP) + aes(fill=RR_cat) + theme_bw() +
scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatio-temporal model") +
theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
# PP plot
rape_PP_st = ggplot() + geom_sf(data=map_st_RR_PP) + aes(fill=PP_cat) + theme_bw() +
scale_fill_viridis(
option = "plasma", name="PP",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP Spatio-temporal model") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
rape_RR_st | rape_PP_st
# Define formula
formula_stInt_dowry = O ~ f(ID_area, model="bym2", graph=adj, hyper=list(
prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
phi=list(prior="pc", param=c(0.5,2/3))
)) +
f(ID_year, model="rw1", hyper=list(prec=list(prior="pc.prec", param=c(0.5/0.31, 0.01)))) +
f(ID_area_year, model="iid", hyper=list(prec=list(prior="pc.prec", param = c(0.5/0.31, 0.01))))
# Fit model
st_int_model_dowry = inla(formula = formula_stInt_dowry, family="poisson", data=data_st_dowry, E=E,
control.compute=list(dic=TRUE, waic=TRUE, config = TRUE, return.marginals.predictor = TRUE),
control.predictor = list(compute=TRUE))
# extract temporal RRs & CIs to be usied for plotting later
RR_stInt_RW_RR = c()
RR_stInt_RW_l = c()
RR_stInt_RW_h = c()
for(i in 1:14) {
RR_stInt_RW_RR[i] = inla.emarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_year[[i]])
RR_stInt_RW_l[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_year[[i]]))
RR_stInt_RW_h[i] = inla.qmarginal(0.975,inla.tmarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_year[[i]]))
}
RR_stInt_RW_dowry = data.frame(RR=RR_stInt_RW_RR, low_CI=RR_stInt_RW_l, up_CI = RR_stInt_RW_h)
# extract spatial RRs
RR_stInt_BYM = c()
for(i in 1:70) {RR_stInt_BYM[i] = inla.emarginal(function(x) exp(x),
st_int_model_dowry$marginals.random$ID_area[[i]])}
RR_stInt_BYM_marg = st_int_model_dowry$marginals.random$ID_area[1:70]
PP_stInt_BYM = lapply(RR_stInt_BYM_marg, function(x) {1-inla.pmarginal(0,x)})
# combine
resRR_PP_stInt = data.frame(RR=RR_stInt_BYM, PP=unlist(PP_stInt_BYM), ID_area=carto_up$ID_area)
# plot (using the process defined above)
breaks = c(min(resRR_PP_stInt$RR), 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(resRR_PP_stInt$RR))
resRR_PP_stInt$RR_cat = cut(resRR_PP_stInt$RR, breaks=breaks, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
resRR_PP_stInt$PP_cat = cut(resRR_PP_stInt$PP, breaks_pp, include.lowest = TRUE)
map_stInt_RR_PP=left_join(carto_up, resRR_PP_stInt, by=c("ID_area"="ID_area"))
# RR plot
dowry_RR_stInt = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=RR_cat) + theme_bw() +
scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatio-temporal model with Interaction I") +
theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
# PP plot
dowry_PP_stInt = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=PP_cat) + theme_bw() +
scale_fill_viridis(
option = "plasma", name="PP",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP Spatio-temporal model with Interaction I") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
dowry_RR_stInt | dowry_PP_stInt
formula_stInt_rape = O ~ f(ID_area, model="bym2", graph=adj, hyper=list(
prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
phi=list(prior="pc", param=c(0.5,2/3))
)) +
f(ID_year, model="rw1", hyper=list(prec=list(prior="pc.prec", param=c(0.5/0.31, 0.01)))) +
f(ID_area_year, model="iid", hyper=list(prec=list(prior="pc.prec", param = c(0.5/0.31, 0.01))))
st_int_model_rape = inla(formula = formula_stInt_rape, family="poisson", data=data_st_rape, E=E,
control.compute=list(dic=TRUE, waic=TRUE))
# temporal RRs & CIs
RR_stInt_RW_RR = c()
RR_stInt_RW_l = c()
RR_stInt_RW_h = c()
for(i in 1:14) {
RR_stInt_RW_RR[i] = inla.emarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_year[[i]])
RR_stInt_RW_l[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_year[[i]]))
RR_stInt_RW_h[i] = inla.qmarginal(0.975,inla.tmarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_year[[i]]))
}
RR_stInt_RW_rape = data.frame(RR=RR_stInt_RW_RR, low_CI=RR_stInt_RW_l, up_CI = RR_stInt_RW_h)
# spatial RRs
RR_stInt_BYM = c()
for(i in 1:70) {RR_stInt_BYM[i] = inla.emarginal(function(x) exp(x),
st_int_model_rape$marginals.random$ID_area[[i]])}
RR_stInt_BYM_marg = st_int_model_rape$marginals.random$ID_area[1:70]
PP_stInt_BYM = lapply(RR_stInt_BYM_marg, function(x) {1-inla.pmarginal(0,x)})
# paste
resRR_PP_stInt = data.frame(RR=RR_stInt_BYM, PP=unlist(PP_stInt_BYM), ID_area=carto_up$ID_area)
breaks = c(min(resRR_PP_stInt$RR), 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(resRR_PP_stInt$RR))
resRR_PP_stInt$RR_cat = cut(resRR_PP_stInt$RR, breaks=breaks, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
resRR_PP_stInt$PP_cat = cut(resRR_PP_stInt$PP, breaks_pp, include.lowest = TRUE)
map_stInt_RR_PP=left_join(carto_up, resRR_PP_stInt, by=c("ID_area"="ID_area"))
c = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=RR_cat) + theme_bw() +
scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatio-temporal model with Interaction I") +
theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
rape_RR_stInt = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=RR_cat) + theme_bw() +
scale_fill_viridis(
option = "plasma", name="PP",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP Spatio-temporal model with Interaction I") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
rape_PP_stInt = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=PP_cat) + theme_bw() +
scale_fill_viridis(
option = "plasma", name="PP",
discrete = T,
direction = -1,
guide = guide_legend(
title.position = 'top',
reverse = T
)) + ggtitle("PP Spatio-temporal model with Interaction I") + theme(text = element_text(size=15),
axis.text.x = element_blank(),
axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))
rape_RR_stInt | rape_PP_stInt
# Define the priors used for sensitivity analysis (to be included in report)
knitr::opts_chunk$set(echo = TRUE)
# Create the table data
prior_settings <- data.frame(
Model = c("Base_PC", "Weaker_Shrinkage", "Stronger_Spatial"),
U = c(0.31, 0.5, 0.3),
Alpha = c(0.01, 0.01, 0.05),
Phi = c("2/3", "0.5", "0.8"),
Description = c(
"Original PC prior",
"Weaker shrinkage: allows more variance, favors less spatial structure",
"Stronger spatial structure: tighter σ, favors more spatial structure"
)
)
# render with borders
sensitivity_priors <- kable(prior_settings, format = "html", escape = FALSE, align = "c",
caption = "Hyperprior Settings for Sensitivity Analysis") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", "bordered"),
full_width = FALSE,
position = "center") %>%
column_spec(1:4, width = "8em") %>%
column_spec(5, width = "30em")
# Function to fit the model with custom PC priors
fit_st_model <- function(data, E, adj, u_val, alpha_val, phi_alpha) {
formula_custom = O ~
f(ID_area, model="bym2", graph=adj,
hyper=list(
prec=list(prior="pc.prec", param=c(0.5/u_val, alpha_val)),
phi=list(prior="pc", param=c(0.5, phi_alpha))
)) +
f(ID_year, model="rw1",
hyper=list(prec=list(prior="pc.prec", param=c(0.5/u_val, alpha_val)))) +
f(ID_area_year, model="iid",
hyper=list(prec=list(prior="pc.prec", param=c(0.5/u_val, alpha_val))))
model = inla(formula_custom, family="poisson", data=data, E=E,
control.compute = list(dic = TRUE, waic = TRUE),
control.predictor = list(compute = TRUE))
return(model)
}
# Define prior settings to test
# Define prior settings to test (updated names)
prior_settings = data.frame(
model_id = c("Base_PC", "Weaker_Shrinkage", "Stronger_Spatial"),
u_val = c(0.31, 0.5, 0.3),
alpha_val = c(0.01, 0.01, 0.05),
phi_alpha = c(2/3, 0.5, 0.8)
)
# Container for results
results = data.frame()
# Loop over models for Dowry and Rape
for(crime in c("dowry", "rape")) {
data_input = if(crime == "dowry") data_st_dowry else data_st_rape
for(i in 1:nrow(prior_settings)) {
cat("Running", crime, "model:", prior_settings$model_id[i], "\n")
model = fit_st_model(
data = data_input,
E = E,
adj = adj,
u_val = prior_settings$u_val[i],
alpha_val = prior_settings$alpha_val[i],
phi_alpha = prior_settings$phi_alpha[i]
)
# Extract hyperparameter summaries
hyp = model$summary.hyperpar
results = rbind(
results,
data.frame(
Crime = crime,
Model = prior_settings$model_id[i],
DIC = model$dic$dic,
WAIC = model$waic$waic,
Prec_ID_area = hyp["Precision for ID_area", "mean"],
Phi_ID_area = hyp["Phi for ID_area", "mean"],
Prec_ID_year = hyp["Precision for ID_year", "mean"],
Prec_ID_area_year = hyp["Precision for ID_area_year", "mean"]
)
)
}
}
## Running dowry model: Base_PC
## Running dowry model: Weaker_Shrinkage
## Running dowry model: Stronger_Spatial
## Running rape model: Base_PC
## Running rape model: Weaker_Shrinkage
## Running rape model: Stronger_Spatial
sensitivity <- kable(results, caption = "Sensitivity Analysis Results for Spatio-temporal Models with Interaction I") %>% kable_styling(bootstrap_options = "striped", full_width = F, position = "center")
# Spatio-temporal model - prepare temporal residual RRs
RR_stRW_dowry$year <- seq(2001, 2014)
RR_stRW_rape$year <- seq(2001, 2014)
RR_stRW_dowry$Crime <- "Dowry Deaths"
RR_stRW_rape$Crime <- "Rape"
RR_stRW_both <- rbind(RR_stRW_dowry, RR_stRW_rape) # Combine
# Spatio-temporal + Int model - prepare temporal residual RRs
RR_stInt_RW_dowry$year <- seq(2001, 2014)
RR_stInt_RW_rape$year <- seq(2001, 2014)
RR_stInt_RW_dowry$Crime <- "Dowry Deaths"
RR_stInt_RW_rape$Crime <- "Rape"
RR_stInt_RW_both <- rbind(RR_stInt_RW_dowry, RR_stInt_RW_rape) # Combine
# Plot - ST model, both crimes
temp_st_both <- ggplot(RR_stRW_both, aes(x = year, y = RR, colour = Crime, fill = Crime)) +
geom_line() +
geom_ribbon(aes(ymin = low_CI, ymax = up_CI), alpha = 0.2, colour = NA) +
ggtitle("Spatio-temporal model") +
theme_bw() +
labs(x = "Year")
# Plot - ST+I model, both crimes
temp_stInt_both = ggplot(RR_stInt_RW_both, aes(x = year, y = RR, colour = Crime, fill = Crime)) +
geom_line() +
geom_ribbon(aes(ymin = low_CI, ymax = up_CI), alpha = 0.2, colour = NA) +
ggtitle("Spatio-temporal model with Interaction I") +
theme_bw() +
labs(x = "Year")
# Combine into one figure
temp_st_stInt <- temp_st_both | temp_stInt_both
temp_st_stInt
For both models
# Prepare data
data_st = left_join(carto_up, data, by = "ID_area")
# Extract interactions
data_st$int_dowry = st_int_model_dowry$summary.random$ID_area_year$mean
data_st$int_dowry_cat = cut(data_st$int_dowry, breaks=c(-1,-0.05,
-0.01, 0.01, 0.05, 1),include.lowest = T)
data_st$int_rape = st_int_model_rape$summary.random$ID_area_year$mean
data_st$int_rape_cat = cut(data_st$int_rape, breaks=c(-1,-0.05,
-0.01, 0.01, 0.05, 1),include.lowest = T)
# Plot interactions for dowry model
dowry_int_map = ggplot() + geom_sf(data = data_st, aes(fill=int_dowry_cat)) + theme_bw() +
scale_fill_brewer(palette = "PuOr") +
guides(fill=guide_legend(title=NULL)) +
theme(text = element_text(size=20),
axis.text.x = element_blank(), axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 5, labeller=labeller(ID_year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005",
"6"="2006","7"="2007","8"="2008","9"="2009","10"="2010",
"11"="2011","12"="2012","13"="2013","14"="2014"))) +
labs(title="Dowry")
# Plot interactions for rape model
rape_int_map = ggplot() + geom_sf(data = data_st, aes(fill=int_rape_cat)) + theme_bw() +
scale_fill_brewer(palette = "PuOr") +
guides(fill=guide_legend(title=NULL)) +
theme(text = element_text(size=20),
axis.text.x = element_blank(), axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 5, labeller=labeller(ID_year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005",
"6"="2006","7"="2007","8"="2008","9"="2009","10"="2010",
"11"="2011","12"="2012","13"="2013","14"="2014"))) +
labs(title="Rape")
# Combine
int_map = rape_int_map / dowry_int_map
int_map
# Exract posterior RRs
data_st$RR_total_emarg_dowry <- mapply(function(area_idx, year_idx, area_year_idx) {
rr_area <- inla.emarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_area[[area_idx]])
rr_year <- inla.emarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_year[[year_idx]])
rr_interaction <- inla.emarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_area_year[[area_year_idx]])
rr_area * rr_year * rr_interaction
},
area_idx = data_st$ID_area,
year_idx = data_st$ID_year,
area_year_idx = data_st$ID_area_year)
# categorise
breaks_total = c(min(data_st$RR_total_emarg_dowry), 0.5, 1, 1.5, 2, 2.5, max(data_st$RR_total_emarg_dowry))
data_st$RR_total_cat_dowry = cut(data_st$RR_total_emarg_dowry, breaks=breaks_total,include.lowest = T)
# Plot for dowry
dowry_stInt_postRR = ggplot() + geom_sf(data = data_st, aes(fill = RR_total_cat_dowry)) + theme_bw() + scale_fill_brewer(palette = "PuOr") +
guides(fill=guide_legend(title=NULL)) +
theme(text = element_text(size=20),
axis.text.x = element_blank(), axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 5, labeller=labeller(ID_year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005",
"6"="2006","7"="2007","8"="2008","9"="2009","10"="2010",
"11"="2011","12"="2012","13"="2013","14"="2014"))) +
labs(title = "Dowry Crime")
data_st$RR_total_emarg_rape <- mapply(function(area_idx, year_idx, area_year_idx) {
rr_area <- inla.emarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_area[[area_idx]])
rr_year <- inla.emarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_year[[year_idx]])
rr_interaction <- inla.emarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_area_year[[area_year_idx]])
rr_area * rr_year * rr_interaction
},
area_idx = data_st$ID_area,
year_idx = data_st$ID_year,
area_year_idx = data_st$ID_area_year)
breaks_total = c(min(data_st$RR_total_emarg_rape), 0.5, 1, 1.5, 2, 2.5, max(data_st$RR_total_emarg_rape))
data_st$RR_total_cat_rape = cut(data_st$RR_total_emarg_rape, breaks=breaks_total,include.lowest = T)
rape_stInt_postRR = ggplot() + geom_sf(data = data_st, aes(fill = RR_total_cat_rape)) + theme_bw() + scale_fill_brewer(palette = "PuOr") +
guides(fill=guide_legend(title=NULL)) +
theme(text = element_text(size=20),
axis.text.x = element_blank(), axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 5, labeller=labeller(ID_year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005",
"6"="2006","7"="2007","8"="2008","9"="2009","10"="2010",
"11"="2011","12"="2012","13"="2013","14"="2014"))) +
labs(title = "Rape")
# Combine posterior RRs & plot
meanRR_post = rape_stInt_postRR / dowry_stInt_postRR + plot_annotation(title="Posterior Mean Relative Risks", theme=theme(plot.title=element_text(size=24)))
meanRR_post
# descriptive analysis
unsmoothed <- usmoothed_map | unsmoothed_smr_map
region_temp_smrs <- (p_dowry | p_rape) +
plot_layout(guides = "collect") +
plot_annotation(title = "Crude Rates of Rape and Dowry Crime in Uttar Pradesh (2001–2014)")
# supplementary figures of all models' smoothed spatial RR maps
rape_all_m <- (rape_RR_s | rape_PP_s) / (rape_RR_st | rape_PP_st) / (rape_RR_stInt | rape_PP_stInt)
dowry_all_m <- (dowry_RR_s | dowry_PP_s) / (dowry_RR_st | dowry_PP_st) / (dowry_RR_stInt | dowry_PP_stInt)
# ST+I RR & PP plots for rape and dowry crime (model of interest discussed in report)
rape_title <- wrap_elements(full = textGrob("Rape", rot = 90, gp = gpar(fontsize = 15))) # Make titles
dowry_title <- wrap_elements(full = textGrob("Dowry Crime", rot = 90, gp = gpar(fontsize = 15)))
rape_row <- rape_title + (rape_RR_stInt + rape_PP_stInt + plot_layout(ncol = 2)) +plot_layout(widths = c(0.1, 1)) # Compile rows
dowry_row <- dowry_title + (dowry_RR_stInt + dowry_PP_stInt + plot_layout(ncol = 2)) +plot_layout(widths = c(0.1, 1))
stInt_m <- rape_row / dowry_row # compile figures
WAIC_df = data.frame(model = c("Spatio-temporal", "Spatio-Temporal + Int"),
WAIC_dowry = round(c(st_dowry_model$waic$waic, st_int_model_dowry$waic$waic)),
WAIC_rape = round(c(st_rape_model$waic$waic, st_int_model_rape$waic$waic)))
colnames(WAIC_df) = c("Model", "WAIC (Dowry)", "WAIC (Rape)")
waic_table = knitr::kable(WAIC_df, caption = "Table 1: WAIC for Different Models and Crimes") %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center")
waic_table
| Model | WAIC (Dowry) | WAIC (Rape) |
|---|---|---|
| Spatio-temporal | 6654 | 7781 |
| Spatio-Temporal + Int | 6318 | 6240 |
# Combine hyperparameters into one table
rape_s_hyper = s_rape_model$summary.hyperpar
dowry_s_hyper = s_dowry_model$summary.hyperpar
rape_st_hyper = st_rape_model$summary.hyperpar
dowry_st_hyper = st_dowry_model$summary.hyperpar
rape_stInt_hyper = st_int_model_rape$summary.hyperpar
dowry_stInt_hyper = st_int_model_dowry$summary.hyperpar
# tidy up into dataframes: spatial
rape_s_df <- as.data.frame(rape_s_hyper) %>%
mutate(Hyperparameter = rownames(rape_s_hyper),Crime = "Rape", Model = "Spatial")
dowry_s_df <- as.data.frame(dowry_s_hyper) %>%
mutate(Hyperparameter = rownames(dowry_s_hyper),Crime = "Dowry",Model = "Spatial")
# spatio-temporal
rape_st_df <- as.data.frame(rape_st_hyper) %>%
mutate(Hyperparameter = rownames(rape_st_hyper),Crime = "Rape",Model = "Spatio-temporal")
dowry_st_df <- as.data.frame(dowry_st_hyper) %>%
mutate(Hyperparameter = rownames(dowry_st_hyper),Crime = "Dowry",Model = "Spatio-temporal")
# ST + interaction
rape_stInt_df <- as.data.frame(rape_stInt_hyper) %>%
mutate(Hyperparameter = rownames(rape_stInt_hyper),Crime = "Rape",Model = "Spatio-temporal Interaction")
dowry_stInt_df <- as.data.frame(dowry_stInt_hyper) %>%
mutate(Hyperparameter = rownames(dowry_stInt_hyper),Crime = "Dowry",Model = "Spatio-temporal Interaction")
# Combine dataframes
hyper_df <- bind_rows(rape_s_df, dowry_s_df,rape_st_df, dowry_st_df,rape_stInt_df, dowry_stInt_df)
# Plot all hyperparameters from all models
hyper_all = ggplot(hyper_df, aes(x = Model, y = mean, color = Crime, group = Crime)) +
geom_point(position = position_dodge(width = 0.5), size = 3) +
geom_errorbar(aes(ymin = `0.025quant`, ymax = `0.975quant`),
width = 0.2, position = position_dodge(width = 0.5)) +
facet_wrap(~ Hyperparameter, scales = "free_y") +
labs(title = "Hyperparameter Estimates by Model and Crime Type",
y = "Estimate", x = "Model") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
hyper_all
rape_hyper <- as.data.frame(st_int_model_rape$summary.hyperpar)
dowry_hyper <- as.data.frame(st_int_model_dowry$summary.hyperpar)
# 2. Clean / Rename Columns for Clarity
# Typically, you have columns like "mean", "sd", "0.025quant", "0.975quant", etc.
rape_hyper <- rape_hyper %>%
mutate(Parameter = rownames(rape_hyper),
CrimeType = "Rape") %>%
select(Parameter, mean, `0.025quant`, `0.975quant`, CrimeType)
dowry_hyper <- dowry_hyper %>%
mutate(Parameter = rownames(dowry_hyper),
CrimeType = "Dowry Deaths") %>%
select(Parameter, mean, `0.025quant`, `0.975quant`, CrimeType)
# Combine into a single data frame
hyper_combined <- bind_rows(rape_hyper, dowry_hyper)
# 3. Create a Nice Table
# This example shows a side-by-side listing of parameter, mean, and credible intervals
# for each crime type.
hyper_table <- hyper_combined %>%
mutate(
# Format the credible interval as a single string
CrI = paste0("(", round(`0.025quant`, 3), ", ", round(`0.975quant`, 3), ")"),
Mean = round(mean, 3)
) %>%
select(CrimeType, Parameter, Mean, CrI)
# 1. Pivot to wide format
hyper_wide <- hyper_combined %>%
pivot_wider(
id_cols = Parameter,
names_from = CrimeType, # Rape vs. Dowry Deaths
values_from = c(mean, `0.025quant`, `0.975quant`),
names_glue = "{CrimeType}_{.value}" # to create columns like Rape_mean, Dowry Deaths_mean, etc.
)
# 2. Create columns that combine the mean and 95% CrI into one string
hyper_wide <- hyper_wide %>%
mutate(
Rape_Est = paste0(
round(Rape_mean, 3),
" (", round(Rape_0.025quant, 3), ", ", round(Rape_0.975quant, 3), ")"
),
Dowry_Est = paste0(
round(`Dowry Deaths_mean`, 3),
" (", round(`Dowry Deaths_0.025quant`, 3), ", ", round(`Dowry Deaths_0.975quant`, 3), ")"
)
) %>%
# 3. Select and rename columns for clarity
select(
Parameter,
`Rape: Mean (95% CrI)` = Rape_Est,
`Dowry Deaths: Mean (95% CrI)` = Dowry_Est
)
# 4. Create a publication-ready table
stInt_hyper = kable(hyper_wide, "html", caption = "Table 2:Posterior Summaries of Hyperparameters") %>%
kable_styling(full_width = FALSE)
stInt_hyper
| Parameter | Rape: Mean (95% CrI) | Dowry Deaths: Mean (95% CrI) |
|---|---|---|
| Precision for ID_area | 5.911 (3.857, 8.543) | 11.343 (7.555, 16.158) |
| Phi for ID_area | 0.849 (0.567, 0.983) | 0.894 (0.675, 0.99) |
| Precision for ID_year | 15.095 (5.969, 30.791) | 52.506 (19.597, 110.876) |
| Precision for ID_area_year | 12.55 (10.765, 14.554) | 43.712 (34.537, 54.93) |
This is done so that figures can be subsequently loaded into the report file and visualised
save(unsmoothed, region_temp_smrs, rape_all_m, dowry_all_m, rape_title, dowry_title, stInt_m, temp_st_stInt, int_map, meanRR_post, waic_table, stInt_hyper, hyper_all, reg_table, sensitivity_priors, sensitivity, file="Figures.RData")